Randomly Evolving Idiotypic Networks: 
Modular Mean Field Theory 
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We develop a modular mean field theory for a minimalistic model of the idiotypic network. The 
model comprises the random influx of new idiotypes and a deterministic selection. It describes the 
evolution of the idiotypic network towards complex modular architectures, the building principles 
of which are known. The nodes of the network can be classified into groups of nodes, the modules, 
which share statistical properties. Each node experiences only the mean influence of the groups to 
which it is linked. Given the size of the groups and linking between them the statistical properties 
such as mean occupation, mean life time, and mean number of occupied neighbors are calculated 
for a variety of patterns and compared with simulations. For a pattern which consists of pairs of 
occupied nodes correlations are taken into account. 

PACS numbers: 87.18.-h, 87.18. Vf, 87.23. Kg, 87.85.Xd, 64.60.aq, 05.10.-a, 02.70.Rr 



I. INTRODUCTION 

B-Lymphocytes express on their surface Y-shaped re- 
ceptors, called antibodies, with highly specific binding 
sites. All antibodies of a given B-cell are of the same 
type, the idiotype. If the antibodies are cross-linked by 
complementary structures, situated e.g. on an antigen, 
the B-cell is stimulated to proliferate and, after a few cell 
cycles, differentiate into plasma cells and memory cells. 
Plasma cells secrete large amounts of antibodies, which 
attach to the antigen and mark it for further processing. 
Useful clones survive, while others, lacking stimulation, 
die [J. 

Complementary structures can be found also on B- 
lymphocytes. B-cells of complementary idiotype may 
stimulate each other, thus the B-lymphocyte system 
forms a functional network, the idiotypic network 
A history and thorough discussion of immunological 
paradigms can be found in [3], cf. also [J]. For reviews 
on idiotypic networks with emphasis on modeling see @ , 
and with focus on new immunological and clinical find- 
ings see @. 

The idiotypic network is an attractive concept for sys- 
tem biologists, but due to their complexity also a chal- 
lenge for theoretical physicists. The size of the poten- 
tial idiotypic repertoire of humans is estimated to exceed 
10 10 Q, the expressed repertoire is of order 10 s Q. In- 
teractions between B-cells of complementary idiotype are 
genuinely nonlinear. Thus, modeling idiotypic networks 
is an inviting playground for statistical physics, nonlinear 
dynamics, and complex systems. Networks, especially 
random and randomly evolving networks, with applica- 
tions in a plethora of different, multidisciplinary fields 
experience great interest in the community of sta- 
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tistical physicists. Computer scientists try to mimick the 
immune system to fight against foreign invaders [l4j ]. 

A minimalistic model of the idiotypic network was pro- 
posed in [l5j where the nodes represent B-lymphocytes 
and antibodies of a given idiotype. The idiotype is char- 
acterized by a bitstring. Populations with complemen- 
tary idiotypes, allowing for a few mismatches, can inter- 
act. In the model, an idiotype population may be present 
or absent. For survival it needs stimulation by sufficiently 
many complementary idiotypes, but becomes extinct if 
too many complementary idiotypes are present. This re- 
flects the log-bell shaped dose-response curve character- 
istic for B-lymphocytes Q. The dynamics is driven by 
the influx of new idiotypes generated by mutations in the 
bone marrow. 

The potential idiotypic network consists of all idiotypes 
an organism is able to generate. Each idiotype v is la- 
belled by a bitstring of length d: 6<jbd-i • • • &1, which is 
the binary address of the node in the network. Two nodes 
v and w are linked if their bitstrings are complementary. 
We allow for m mismatches, i.e. their Hamming distance 
must obey dn(v, w) ^ d—m. These nodes and links build 
an undirected graph , the base graph. Each node has 

the same number of neighbors, k = X^fcLo it) 1 ^ ne ex ~ 
pressed idiotypic network is only a part of the potential 
network. 

New idiotypes generated in the bone marrow are in- 
troduced by occupying empty nodes randomly. Occu- 
pied nodes are selected to survive if they receive suffi- 
cient stimulus, i.e. if the number of occupied neighbors 
is within an allowed window. To be specific, the rules for 
(parallel) update are 

(i) Occupy empty nodes with probability p 

(ii) Count the number of occupied neighbors n(dv) of 
node v. If n{dv) is outside the window [tL,tu] > set 
the node v empty 

(iii) Iterate. 
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The model has a minimal number of parameters, 
namely the length of the bitstring, the allowed number of 
mismatches, upper and lower thresholds of the window, 
and the influx probability of new idiotypes. 

We can consider our model system as a descendant of 
Conway's game of life 16|. There, on an infinite regular 
2d lattice the sites can take value or 1. The system 
is updated in parallel in discrete time. If the number 
of living neighbors lies between a lower and an upper 
threshold, a site becomes populated or survives in the 
next step. The dynamics is entirely determined by the 
initial configuration. There is continuous interest in the 
highly complex properties of game of life, for a recent 
status report see [171 ] . 

Schulman and Seiden fl8j investigated a probabilistic 
version of the game of life, where the update rule is mod- 
ified in two ways. Sites for which a window rule is ful- 
filled will be occupied or will survive with given probabil- 
ities. Furthermore, a site can be occupied or survives in 
a stochastic way, parametrized by a temperature T such 
that for T = the modified window rule applies but plays 
no role for T — > oo. Starting from random initial condi- 
tions, simulations show a sharp transition of the global 
mean occupation when T is increased. The high tem- 
perature phase is well described by a mean field theory 
which fails however to reproduce low temperature results. 
Excluding sites without occupied neighbors improves the 
agreement of theory and simulation for low temperatures. 
A second order mean field theory was proposed in (l9j |. 

Gutowitz et al. [2(| developed a local structure theory 
for cellular automata on regular lattices, which considers 
the evolution of joint probabilities of sets of neighbor- 
ing cells (block configurations) and thus include corre- 
lations. Application to the problem of Schulman and 
Seiden yields results which agree with their simulations 

a 

Bidaux et al. [22| considered a binary probabilistic 
cellular automaton with a totalistic update rule involv- 
ing 8 neighbors on regular lattices in d = 1,...,4. If 
the window rule is fulfilled, it applies with probability p. 
Simulations with random initial conditions show a tran- 
sition of the global mean occupation at a critical value 
p c , which is continuous for d = 1 and first order for d > 1. 
A mean field theory describes qualitatively a first order 
transition. An overview of mean field theories of cellular 
automata including the probabilistic game of life can be 
found in 23] 

In the context of network models there are numerous 
mean field approaches aimed to describe degree distribu- 
tion, clustering coefficient, average shortest distance be- 
tween two nodes, mean populations, e.g. in susceptiblc- 
infectious-recovered models of epidemic spread, and 
other characteristics. Many references can be found in 
the recent monograph by Newman [l3| , cf . also [TJ H3~ 

m. 

Gleeson and Cahalane [27j investigate cascades of acti- 
vation in a model with a threshold dynamics on a Poisson 
random graph with average degree z. Starting with few 



active nodes, neighbors become permanently active if the 
fraction of active neighbors exceeds a threshold, drawn 
for each node from a given distribution. The fraction of 
activated nodes in the nth update step is determined re- 
cursively. The final fraction of active nodes p is related 
to the fixed points of this recursion. For Gaussian dis- 
tributed thresholds with mean R this fraction as a func- 
tion of z may change in a continuous or discontinuous 
way depending on R. 

In this paper we deal with a network of complex archi- 
tecture which emerges as a result of a random evolution. 
The architecture is build of groups of nodes, the mod- 
ules, which share statistical properties. To describe this 
architecture, a global mean field theory is obviously not 
appropriate. Depending on the parameters, different ar- 
chitectures can occur. In a previous paper [281 ] we have 
reached a detailed understanding of the building princi- 
ples of the architectures, which allows to calculate the 
number and size of the groups and their linking. Here 
we develop a modular mean field theory to calculate the 
statistical properties for given architectures. 

In the next section we shortly sketch the building prin- 
ciples [HI and collect the results needed to develop the 
theory. In Sec. IIIII we determine the evolution equation 
for the mean population of the groups for a given ar- 
chitecture, which properly takes into account the update 
rule of the model, random influx and the window rule. 
The fixed point of this evolution equation gives the sta- 
tionary mean populations and allows to compute also the 
mean life time and the mean occupation of the neighbor- 
hood. In Sec. IIVI we consider two groups of nodes for 
which the mean field theory considerably simplifies: Sin- 
gletons, which are essentially isolated nodes, and groups 
of self-coupled nodes, core groups, in static patterns. In 
Sec. |V] we extend the mean field theory by including cor- 
relations, which is necessary to describe 2-cluster pat- 
terns. The results of the mean field theory for the general 
case are compared with simulations for a range of the in- 
flux parameter p in Sec. IVII Some probabilistic aspects 
of the stability of patterns are discussed in the appendix. 



II. ARCHITECTURE OF PATTERNS 

Simulations of the model for one and two allowed mis- 
matches revealed that the system evolves for typical pa- 
rameters towards a complex functional architecture |15| . 
Groups of nodes were identified which share statistical 
properties such as the mean occupation, the mean life 
time and the mean number of occupied neighbors. Also 
the size of the groups and the linking between them have 
been determined. With increasing influx of new idio- 
types transitions between architectures of different com- 
plexity are observed. For small influx static patterns are 
found. For an intermediate range of the influx a station- 
ary dynamic architecture is observed which is the most 
interesting one. It includes a densely connected core, 
a periphery, and isolated nodes (singletons), resembling 
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the notion of central and peripheral part of the biolog- 
ical network [2!| [3(| ■ For larger influx the architecture 
becomes more irregular. 

In [2^, [3l[ an analytic description of the general build- 
ing principles of these architectures was proposed. It 
allows to calculate number and size of groups and their 
linking for a given architecture. Ideal static patterns, i.e. 
patterns without defects which persist without influx, are 
completely characterized. 

For a given architecture, i.e. a pattern, the nodes can 
be classified according to the values of bits in determinant 
positions common to all nodes. Different patterns are 
characterized by the number du of such positions. The 
entries in these determinant positions decide to which 
group a node belongs. The pattern can be built by reg- 
ular arrangement of elementary building blocks, which 
are hypercubes of dimension dm- We call these build- 
ing blocks pattern modules. The concept is explained in 
detail and many examples are given in [2&| . 

For a given pattern, i.e. a given c?m, the number of 
groups and their sizes can be calculated combinatorially 
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FIG. 1. (Color online) A visualization of the 12-group 
architecture, cLm = 11, observed in simulations for p ~ 
0.025 ... 0.045 on a base graph G@ with [t L ,tu] = [1,10] 
in agreement with Eqs. {TJ and ([2]). The size of the boxes 
corresponds to the group size. The lines show possible links 
between nodes of the groups. 



and the number of links L g i of a given node v g G S g 
nodes in Si is 



to 



X 



d-d M \ ( 9-1 
j J Vi+max(0, -A gl ) 

d-M -g + l 
i + max(0, A g i) 



l(j + 2i+\A gl \ <m) (2) 



where A g i = djw— 9~ i+2. 

For example, on a base graph G\ 2 ' with [£l,%] = 
[1, 10] we observe an architecture of three groups. It has 
two determinant bits and is described by a pattern mod- 
ule of dimension dju = 2. Persisting occupied nodes in 
the first group form 2-clusters, potential hubs, the second 
group, occasionally link together several 2-clusters, and 
the nodes in the third group are stable holes. The links 
are given by the corresponding matrix L = (L g i). 

On the same base graph, for intermediate influx we 
also find a dynamical 12-group architecture. Its architec- 
ture, visualized in Fig. [TJ is based on a pattern module 
of dimension dia = H- The corresponding link matrix 
is given in Table HI We find two large core groups with 
links to almost all other groups, two peripheral groups 
connected to the core, groups of stable holes, which sep- 
arate the singletons from the central network. 

In both cases the architecture is stationary, simulation 
and analytic predictions agree perfectly. 

The occupation of nodes is in general not permanent 
but fluctuates due to the influx while the architecture 
persists. For a given pattern the mean occupation of the 
groups varies systematically with the influx p. The mean 
occupation and other statistical quantities can be com- 
puted in a modular mean field theory. We consider the 



TABLE I. Linkmatrix and qualitative classification of groups 
for the 12-group architecture of Fig. [T] Missing entries are 
zero. 
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nodes as situated in a mean field exerted by the neigh- 
boring groups characteristic for a given architecture. The 
required link matrix L in our model can be calculated 
or obtained in simulations. Other modular architectures 
could be treated, too. 
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III. MEAN FIELD APPROACH 

A. Mean Occupation 

In a mean field approach we assume that all nodes are 
independently occupied with a probability Prob(n(i> g ) = 
1) , v g S S g , characteristic for the group the node belongs 
to. Since a node is either occupied or unoccupied, n(v g ) £ 
{0, 1}, we have 

( n ( v g)) = n(vg)Prob(n{vg)) 

n(v g )=0,l 

= Prob(n(v g ) = l) = n g . (3) 

We consider the evolution of a given set of mean occupa- 
tions n = (rii, . . . , nd M +i) to a new set n' induced by the 
update algorithm 

n' = f(n) , (4) 

where / = (fi(n), . . . , fd M +i{n)) . The fixed point n* 
of Eq. @ solves the self-consistency condition that we 
require for the stationary state of our model system 

n* = f(n*) . (5) 

Having obtained the fixed points n*, we can compute 
the mean life time, see Eq. (j2"8")l below, and the number 
of occupied neighbors by 

d M +i 

(n{dv g )) = Y L gi n i ■ ( 6 ) 
1=1 



B. Update map 

1. General approach 

The update map f(n) is constructed following the 
steps of the update algorithm, random influx and appli- 
cation of the deterministic window rule. To be specific, 
we consider a node v g s S g , with g = 1,... , + 1 , and 
its neighbourhood dv g . Then the update algorithm says 
that an empty node v g is occupied by the influx with 
probability p. Next, the window condition is applied: 
An occupied node v g survives if the neighbourhood ful- 
fills the window condition, otherwise it is emptied. We 
denote the occupation of a node v g after the influx by 
n(vg) and the occupation of v g after the update by n(v' g ) . 
A microscopic configuration of dv g after the influx is de- 
noted by C — C(dvg) = («(w)) ligai , which is a list of 
occupations of all neighbors of v g . Then we have 

Prob(n(v' g )) = Y T updatc (n(v'g)\n(vg))PTob(n(vg)) , 

n(v g )=0,l 

(7) 



where 

T up date(n'|n) = (<*n',i<W + pS n ' ,lS n fl)W(dv g ) 

+5 n >, 5 nA (l-W(dv g )) (8) 

is the transition matrix of the update of n(v g ) to n(v' g ) 
depending on the window condition 

W{dv g ) = t(t L < n(dv g ) < t v ) • (9) 

Specifying Eq. © to n(v') — 1 we obtain from Eqs. ([3]), 
0, and © 

n'g = [n g +p(l-n g )]W(dvg). (10) 

Averaging over all configurations {C} leads to 

n; = K+p(l-n s )]P g w , (11) 

where P g w = Prob(W(9w g ) = 1) is the probability that 
after the influx the neighbours of v g fulfill the window 
condition, 

P g w = (W(dv g )) {d} = £ W(9« 9 )Prob(C) 

{C} 

= W(d5 9 ) Y T inRux (C\C)Pmb(C) . (12) 
{C} {c} 

In the next subsections IIII B 21 and IIII B 31 we determine 
the probability Prob(C) of a microstate C = C(dv g ) and 
the transition probability Ti n flux(C|C) for the transition 
from C to C induced by the influx. In subsection IIII B 41 
Pj^ is explicitly determined. 

2. Probability of a microstate 

We introduce the notation (dv g )i = dv g DSi for the set 
of neighbors of v g belonging to group Si, \(dv g )i\ = L g t. 
A microscopic configuration of {dv g )i is denoted by C; = 
{n(w)) we (g v v . Figure[5]shows an example configuration. 
Furthermore, a microscopic configuration C\ with fc; = 
n((dv g )i) occupied nodes is denoted by Ci\k r Such a 
configuration has probability 

Prob(Q| fcl ) = (n ; ) fei (l-n ; ) L9l - fcl . (13) 

There are ( fc 9i ) equivalent microconfigurations with the 
same number of occupied nodes fc;. Multiplying this 
number yields, of course, the binomial distribution. The 
multiplicity has to be taken into account when calculat- 
ing the average occupation of (dv g )i 



{n((dv g )i)) {Cllki } = 




Prob(C;| fc ,) = L g im . 



(14) 

The probability of a microscopic configuration Ct^x = 
Uz Ci\k t of the whole neighbourhood dv g is determined 
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all mesostates C^} and C{fc;} with occupation numbers 

{ki} and {ki} we have to account for the multiplicities 
derived above. This leads to 



E 



ki=0 



d M +l 



1 (t L < n{dv g ) s$ it;) 



1=1 



dVn 



n|(tt)a: 



xT influx (C Z | fei Qi/cJProbCQifcJ . 



(18) 



FIG. 2. A diagram of a microscopic configuration of v g 's 
neighborhood. Filled and empty circles represent occupied 
nodes and unoccupied nodes, respectively. 



by the occupation numbers {ki} , I = 1,...,c£m + 1- It 
factorises as 

d-M+ 1 

Prob(C {M )= [] Prob(Q|0- (15) 
i=i 



3. Random influx 



The probability of the transition between two mi- 
crostates C;| fe; and Cn^i IS 



where n(dv g ) = X)f=i +1 ^- The factors 

Tinfl ux ((5ii| fel \Ci\ kl ) and Prob(C;| fci ) are given by Eqs. (TTB]) 
and US]). Since by definition (~ 9 'J* 1 ) = for h~h < 
and ki — ki > _L g ; — fc; , it is not necessary to explicitly 
write the factor l(fc; — ki ^ 0), in contrast to Eq. (|16[) . 



We obtain 



nW _ 

a ~ 



E 



d.M+1 k, 

><nE 

Z=l fe,=o 



dM+l N 

\ i=i / 



Lgl - h 

h - h 



Lgl 



(19) 



i'o/ — ki 



T'influx(C;|fc i |C;|fc z 



fci— ki 



(l-p) L s i- k il(k l -ki^0) 



(16) 

Here, /c; — ki is the number of empty nodes in CW which 
become occupied, L g i — fc; is the number of nodes of 
CW which remain empty. The factor l(fc; — k\ ^ 0) 
is introduced the number of occupied nodes can not de- 
crease during the influx. There are ( L ~ l , k ') different mi- 

v ki—ki ' 

croconfigurations which can be reached from one 

microconfiguration CW with the same probability, cf. 

Eq. (fTBT) . Again, there are ( L k 9 ') equivalent microconfig- 
urations C;|fc, • Considering the whole neighbourhood we 
have 

djvf+l 

?l„flux(C{A; i }|C{ fei }) = T in fl ux (C;| fci \Ci\k,) ■ (17) 

1=1 



4- Window rule 

After the random influx which leads to a transition 
from C = C(dv g ) to C = C(dv g ) it is tested whether 
or not the window condition Eq. (JUJ) W(dv g ) = 1(£l < 
n(dvg) ^tjj) is fulfilled. The probability that the window 
condition is fulfilled is given by Eq. (fT2")l . Replacing the 
sums over all microstates of C and C by the sums over 



which holds for all groups g = 1 , . . . , <1m + 1 • We can 
simplify Eq. (|19[) observing 



Lgl 

h 



Lgl - h 

h - h 



L g i\ (ki 
k. 



(20) 



Applying the binomial formula 



E 

ki=0 



' ~. \ki—kikj—ki 

n, (1 — ni) p 



[ni + (1 - ni)p] kl , 
(21) 



we can carry out the sum over the k\ in Eq. (|19[) and 
obtain 



E 

fc;=0. 
d M +l 

n 

i=i 



d M +i 



dM+l 



1=1 

L 



i=i 



~) {ni + (l-ni)p] k < 



:[l-ni- (1 -m)p] Lg 



Lol-ki 



(22) 



Since the mean occupation of some node in Si after the 
influx is ni + (1 — ni)p = ni , Eq. (|22[) can be written in 
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a compact way as 



E 

k,=0 



d M +l 



d M +l 



1 = 1 



1 = 1 



x II (j-})xf(l-m) L ° l - kl , (23) 

which has an obvious intuitive meaning: The probabil- 
ity that ki of the L g i nodes of (dv g )i are occupied is 
( ~ 9l )nf ! (l - ni) L ^'- li ' and P g w is obtained by multiply- 
ing over all groups and summing over all occupations 
obeying the window condition. Thus, all quantities in 
Eq. (jlip are determined and we are able to calculate the 
mean occupation of all groups ra, e.g. by iteration of 
Eq. g|). 

C. Mean life time 

In a similar way, we are now able to calculate the mean 
life time (expectation of life) of an occupied node v g in 
group S g . Its expectation value is defined as 



(r(v g )) = r g = o-Prob(r(v ff ) = a) 



(24) 



cr=0 



where Prob(Y(i; g ) = a) is the probability that an oc- 
cupied node v g € S g remains occupied in a subsequent 
steps and disappears in the following step, 

Prob(r(u 9 ) = a) = (25) 
Prob(n CT+ i(u g ) = 0, n„(v g ) = l, . . . , ni(v g ) = l \ n (v g ) = l) 

This can be expressed in terms of the update transition 
matrix T up dnte(n' \n) introduced in Eq. (|8]l. We use the 
shorthand T( CT )(1|1) = T updatc (n a (v g ) = 1 | n„-i(v g ) = 1) 
and T( CT )(0|1) defined accordingly and write 

T( CT+1 )(0|1)T^(1|1) . . .TW(1|1) (26) 

= (i - w(a^ CT+1) ))w(<9^' T) ) . . . w(dv g V) , 

where dv^ denotes the neighbourhood of v g after the 
influx in the er-th step of the iteration. 

Now we take the average of Eq. (|2"B"|) over all possible 
configurations . Assuming that the configurations in 
consecutive steps are independent, the average of Eq. ([26 
factorizes. In the steady state (' 
independent of the time step a and given by Eq 
This yields 



'^9 )){c(")} — Pg* * s 



Pvob(r(v g )=a) = (1 - R W )(P 9 W ) C 



(27) 



Observing £~o a(P g w 
obtain 



p W/(1 _ pW)2 we finally 



(28) 



IV. SIMPLE SPECIAL CASES 

In this section we discuss the numerical calculation of 
statistical network properties and consider special cases 
which allow essential simplifications. The general case is 
considered in Sec. IVI1 

We consider the model on the basegraph with 
parameters [t^ , tu] and p. We formulate the mean field 
theory for an architecture with module dimension d,M- 
The link matrix L = (L g i) is given by Eq. ([2]). The mean 
occupations of the groups in the steady state are the fixed 
points n* of Eq. ([5|). In general, Eq. ([U is a system of 
c?m + 1 multivariate polynomial equations in the form of 
Eq. (fTTj) . The factor is a polynomial in p of maximal 

order k, where k = Y^Lq ( •) is the degree of a node on 
the base graph. P g is also a multivariate polynomial in 
n;, I = 0, . . . , (Lm + 1- The maximal exponent of each 
ni is L g i , such that the order of the multivariate polyno- 
mial is k, because L g i = n. Note, that it is not 
possible to neglect higher orders of the polynomial, since 
the binomial coefficients may be very large. That is, in 
general we have to solve very complex equations. 

In the following we treat groups of nodes for which 
Eq. ([5]) can be simplified in good approximation exploit- 
ing structural properties of the link matrix and properties 
of the groups of singletons and stable holes. 



A. Singletons 

Singletons are surrounded only by stable holes. They 
occur in most patterns, in static patterns as well as in dy- 
namic ones. Therefore, singletons are important. Their 
treatment is simple, since they decouple in very good 
approximation from the rest. 

Typically, stable holes are surrounded by more than tjj 
occupied nodes, such that their occupation is zero after 
each update step. An occupied stable hole is a very rare 
event, and in good approximation we can assume that all 
stable holes are empty. Singletons are surrounded by k 
stable holes. The probability that k out of n empty stable 
holes are occupied by the influx is (£)p fe (l— p) K ~ k . Thus, 
the probability that the window condition is fulfilled is 
given by 



tu 

p w = V 

sing / j 
k=t L 



P k (l- P ) K 



— k 



(29) 



This follows also directly from Eq. (1221) for singletons set- 
ting all hole groups empty. Note, that Pj^ g om y depends 
on p, and PX g ~ K P f° r small p. For singletons Eq. (fTTj) 
has the unique solution 



P P W 

f sir. 



(30) 



The mean occupation n* ing (p) is plotted in Fig. [3] and 
compared with simulation results. 
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0.04 0.06 
Influx p 



(2) 

FIG. 3. Mean field results (line) for singletons on graph G\ 2 
with [ti, tu] = [1,10] compared with simulation (symbols). 
For small p, n* ing = up 2 , for p 2; 0.075 the patterns in the 
simulation become transient, such that no group averages are 
available. Singletons occur in various patterns, cf. Fig. [5] 



B. Occupied Core in Static Patterns 

Self-coupled nodes in static patterns have only stable 
holes and members of their own group as neighbors. They 
appear for instance in cIm — 2, 4, and 6 patterns, cf. 
28]. It is characteristic for these patterns, that the self- 
coupled nodes have a high mean occupation and thus 
suppress the occupation of the stable holes. Supposing a 
high occupation of the self-coupled nodes we can consider 
stable holes as unoccupied. Then, instead of the system 
of Eqs. we have only one independent equation for 
the respective occupied core group. 

As an example we consider the cLm — 4 pattern. S\ are 
singletons, the group of self-coupled nodes is S 2 - The sta- 
ble hole groups S3, S4, and S5 are suppressed by the occu- 
pation of S 2 if n 2 > n!f ltlcal = %/min(L 3! 2,£4,2,£5,2) = 
0.56. With this constraint we can simply consider the 
union S s \l = S3 U S4 U S5 and put n$ = n.4 = n$ = n s h = 
for all groups of stable holes. This simplifies the link ma- 
trix as given in Fig. |4j The nodes in Si, singletons, are 
treated as given above, n\ = n* ing . For the self-coupled 
nodes in 5*2 we have a polynomial of 4th order in n 2 and 
of 80th order in p, 

3 76 



A ={n 2 +p(\-n 2 )) X] !(! < k 2 +k sh < 10) 
3 



k 2 =0 fe Bh =0 

\k 2 



3-fe 2 



I fc j("2 +p(l-n 2 )) fc2 (l -n 2 -p(l-n 2 )) 
7fi \ 

, JpWl-p) 76 -**. (31) 

KshJ 

In Fig. [5] we plotted n' 2 — n 2 over n 2 for representative 
choices of the influx parameter p. For small values of 
p, n 2 is very close to n 2 such that the common way to 
plot n' 2 over n 2 is inconvenient. The roots of n 2 — n 2 = 
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FIG. 4. The link matrix can be simplified if the stable holes 
are practically always empty after the update, i.e. for n2 3> 
^critical ^ gee ^ ex ^^ Then we lump all stable hole groups S3, S4 
and S5 to S s h- 
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0.01 



-0.01 



-0.02 



i> # 




0.2 0.4 



FIG. 5. Fixed points for self-coupled nodes in the d,M — 4 
pattern on graph with [ti,ttr] = [1,10]. The lines give 
the function n' 2 — n2 from the mean field theory for p — 
{solid), p = 0.015 (dashed), p = 0.025 (dotted), and p = 0.035 
(dash-dotted) . The symbols mark stable (•) and unstable (o) 
fixed points. Note, that the fixed points n£ = and n% = 
0.039 are much smaller than ri2 rltlcal , which is a range where 
the equation is actually not valid. 

give the fixed points n 2 of Eq. (f5Tj) . Fixed points with 
|d/dn2n 2 l < 1, i- e - w ith 



2<— (ni-n 2 )<0. 



(32) 



are stable. 

The case of cIm — 6 is treated in the same way. For 
both d,M — 4 and 6 the stable fixed points are in very 
good agreement with the simulations, cf. Fig. [5] 

The case of (1m = 2 can not be treated this way, how- 
ever. It is easy to see that for p = the respective it- 
eration equation for the occupied group Si is n\ = n\, 
which has the fixed points n* = and 1. The fixed point 
n* = 1 corresponds to the perfect 2-cluster pattern but 
it is unstable. Generically, the iteration converges to the 
stable fixed point n* = 0. Also for p > the iteration will 
not produce the correct result of the simulations, even if 
we consider the full set of equations for all groups, cf. 
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0.7 - 

0.01 0.02 0.03 0.04 

Influx p 

FIG. 6. Mean field results for nodes of S2 in a d,M = 4 pat- 
tern (upper curve) and nodes of group S3 in a d m = 6 pattern 
(lower curve), both on graph G^J with [tL,tu] = [1, 10] com- 
pared with the respective values from simulations (symbols). 
Lines end at the value of p where the patterns become unsta- 
ble in the simulations, cf. Fig. [8] (right). 

Eq. ([5]). The reason is simple: in the 2-cluster pattern an 
occupied node has only one occupied neighbor, which ob- 
viously obstructs the mean field approach. Correlations 
between the two occupied neighbors are not negligible, 
they are treated in the next section. 

V. 2-CLUSTER PATTERN WITH 
CORRELATIONS 



dw dv 

FIG. 7. Nodes v and w from Si forming a 2-cluster and their 
disjoint neighborhoods dv and dw. The line width of the links 
corresponds to the correlation strength. In the ideal pattern 
v and w are occupied and the other nodes are empty. 



measured for all i and j for the 2-cluster and the 8-cluster 
pattern, respectively. The correlation in the 2-cluster is 
the strongest and thus indeed relevant. Assuming inde- 
pendence gives qualitatively false results, as explained 
above. 



TABLE II. dj for the 2-cluster pattern, p = 0.0366. The 
corresponding value for G11 calculated in the mean field 
approach below is (vw) — (v) 2 = 2.1 • 10 -2 . 





Si 


s 2 


s 3 


Si 


2.2 ■ 10" 2 


-2.6 ■ 10 -4 


5.4 ■ 10" 20 


s-2 


-2.6 • 10" 4 


-3.1 ■ 10" 6 


1.2 ■ 10~ 19 


S3 


5.4 • 10~ 20 


1.2 ■ 10 -19 


0.0 



A. Mean occupation 

In any 2-cluster pattern the occupied nodes form pairs 
of mutually stimulating nodes, cf. Fig. [7] Their survival 
significantly depends on the presence of the respective 
partner, they are strongly correlated. This holds for all 
2-cluster patterns on base graphs with arbitrary numbers 
of allowed mismatches m. The pattern module of a 2- 
cluster pattern on G ™ is of dimension = m and the 
group of 2-cluster nodes shall be denoted by Si in all 
these cases. 

We measured in simulations the correlation between 
nearest neighbors in different patterns. The correlation 
is quantified by the two-point connected correlation func- 
tion 

G c (v,w) = (n(v)n(w)) — (n(v))(n(w)) . (33) 

If G c (v,w) > 0, v and w are correlated, if G c (v,w) = 0, 
they are uncorrelated, and if G c (v, w) < 0, they are 
anti-correlated. We determined G c (vi,Wj) for V{ 6 5, 
and Wj £ (dvi)j = dvi n Sj a nearest neighbor of U{ 
and an element of Sj. The average over all members 
of Si and their neighbors in Sj is denoted by Gij — 
{{G c (Vi,Wj)) VieSi ) W:j ^ dvi ) r Tables HI and HII] show CI- 



TABLE III. Gij for the 8-cluster pattern, p = 0.0366. Entries 
are missing if there exist no nearest neighbors in these groups. 



Si S2 S3 S4, S5 

51 -4.9-10" 6 1.3-10 -23 1.3-10 -23 

5 2 -7.0-10" 3 -9.2-10- 6 -3.0-10~ 22 -3.0-10~ 22 

5 3 -4.9-10~ 6 -9.2-10~ 6 -5.0-10- 10 0.0 0.0 
£4 1.3-10~ 23 - 3.0-10 - 22 0.0 0.0 

S 5 1.3-10 - 23 - 3.0-10 - 22 0.0 



Consequently, we must not factorize the joint proba- 
bility Prob(n(w), n(w)) for nodes v and w of a 2-cluster. 
In the following we abbreviate n(v) simply by v, and 
the occupation after the influx n(v) and after the up- 
date n(v') by v and v', respectively. In general it holds 
P(v) = P(v, w). For symmetry reasons (v) = (w) 
and also P(v,w) = P(w,v). And, of course, there is the 
normalization ^2 V P(v, w) = 1. 

As a consequence of the strong correlation between 
nodes of a 2-cluster the update rule n' cor = / CO r(«cor) 
is not only a function of n\ = (v) = P(v = l), but also of 
the pair correlation (vw). 

"cor = ((«), (vw},n 2 ,n 3 ) . (34) 
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We now construct the update map / cor . Therefore, we 
once more determine the transition matrix, cf. Eq. ([5]). 
this time for the pair v and w 

P(v', w') = T updatc {v' 7 w'\v, w)P(v, w) . (35) 



It is constructed from the two subsequent steps of influx 
and application of the window rule 

T updatc (v',w'\v,w) = 

^ Twindow(V, w'\v, w)T in fi ux (v, w\v, w) . (36) 

v,w 

The influx step, which is still independent for the nodes 
v and w is defined as 

71nfiux(w, w\v,w) = r influx (u|w)T in fl ux («}|u;) = 

(Sv,i8v,i + ^,i<5ti,oP + #c,o<^,o(l - p)) (37) 
x {5w,i^w,x + <W<VoP + $w,o8w,o(l - P)) ■ 

It is now straightforward to specify P(v, w) = P(C, w) = 

J2v,w T infiux(v,w\v,w)P(v,w) to 

P(0,0)=P(0,0)(l-p) 2 , 

P(0, 1) = P(0, 0) P (1 -p) + P(0, 1)(1 - p) , 

P(l, I) = P(0, 0)p 2 + 2P(0, l)p + P(l, I) , (38) 

where we used P(1,0) = P(0, 1). Similarily, for the win- 
dow rule we have 

T window (v',w'\v,w) = 
(«t,',i«e,iW(05|t5) +6 <0 Sv,i (1-W(dv\w)) +<V,oM 
x (^i^iW^ltflJ + tf^o^iCl-WCSffilS)) 

+ ^tu'.ofe.o) ) (39) 



where 



W(0« | w = n) = 1 ft L < n + ^ fc ; < tu^j 



(40) 



Here, we sum over the groups S% and S3 and the occupa- 
tion of the partner node w in S\ is explicitly taken into 
account. Note also the explicit dependence of the two 
factors in Eq. (I3T))) on v and w. 

Application of the window rule leads to 

P'(0, 0) = P(0, 0) + P(0, l)(l-W(<9w|0)) 
+P(l,0)(l-W(9u|0)) 

+p(i, i)(i-w(0t&|i))(i-w(0e|i)) , 

P'(0, 1) = P(0, l)W(<9w|0) 

+P(I, l)(l-W(5S|l))W(0u)|l) , 
P'(l,0) = P(l,0)W(5f5|0) 

+P(I, l)W(9u|l)(l-W(9w|l)) , 
P'(l, I) = P(l, l)W(0S|l)W(0u;|l) . (41) 



Averaging Eq. (|40l) over the neighborhood microconfigu- 
rations C(dv\w) of the 2-cluster nodes yields 



(W(dv\n)) {c{a ^)} = Pr 



(42) 



the conditional probability that the window condition is 
fulfilled for v given that the partner node w has occupa- 
tion n. Explicitly, in analogy to Eq. ([23]) 



pW _ 
Mln — 



1 f,<Il 



E 

fei=0 

rii 

1=1 



1=2 

Lxi 

h 



ni) 



1=1 



\L 1 ,-k, 



(43) 



Of course, it holds that 
J2 P^vob(w=n) = P^ (l-H) + P^(w) = P, w . 



n=0,l 

The four pair probabilities P(v,w) are not indepen- 
dent. They obey the symmetry condition and are nor- 
malized. We can express them by (v) and (vw) as 

P(1, 1) = vwP(v, w) — (vw) , 

P(1,0) = (v(l-w)) = (v) - (vw) , 
P(0, 1) = {(l-v)w) = (w) - (vw) , 
P(0,0) = ((l-v)(l-w)) = l-(v)-(w) + (vw), (45) 

where we used (v) = P(v = l) = P(l, I) + P(l, 0). Thus, 
we can use (v) and (vw) as independent variables. The 
update rules for them obtained from Eqs. (|38|) . (14T1) . and 
(1121 are 

(«>' = [-(« w )(i-p) 2 + (v)(i-p)(i-2 P )+p(i-p)] ?; 

+ [(^)(l-p) 2 + (v)2p(l-p) +p 2 ] P* , (46) 

(vw)' = [(^)(l-p) 2 + («>2p(l-p) + P 2 ] (P^) 2 . (47) 

Neglecting the pair correlation reproduces the simple 
mean field equation, of course. 

Finally, the mean occupation of all other groups in 
the 2-cluster pattern are obtained from Eq. (fTTj) using 
ni = (v) from Eq. (gBJ. 

B. Mean life time 

In order to calculate Prob(r(ui) = a), cf. Eq. 
we need the transition probablilities T^(v'\v) 
P(v', v)jP(v) again. However, for the case of strongly 
correlated nodes we have to consider the nodes in a 2- 
cluster v and w as a pair. 

P(v',v) = P(v',w';v,w) , (48) 

w,w' 

P(v', w'; v, w) = T updatc (v', w'\v, w)P(v, w) . (49) 
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Thus, the transition probability corresponding to 
T^iv'lv) in Eq. is now given by 

T cot ( v '\ v ) = E T updatc (vt,Wt\v t ^ 1 ,wt-i)P(v,w)/P(v) , 

W,W f 

(50) 

where P(v) = P(v, w). 

Equation ([50)) can be applied to formula (|26|) in order 
to obtain Prob(r(ui) = a). Explicitly we have 

T c ( or 3 (l|l) = {[(l-p)W(9u|0)+pW(9w|l)]P(l,0) 

+W(dv\l)P(l,l)}/P{l), (51) 
T C W(0|1) = 1-T C W(1|1). (52) 

Using that the two neighbourhoods dv\w and dw\v 
are disjoint and assuming that their configurations after 
the influx (dv\w) in subsequent time steps a are in- 
dependent, we average over all these configurations. This 
replaces all W(dv\n) by P^ n given by Eq. (|4"3")l . Express- 
ing the P{v, w) by the independent variables (v) and (vw) 
we obtain 

<T cor (l|l)) {cM(a{A „ )} = [(l-p)P^ + P P^} (53) 
+ l(l-p)P™+(l-p)P™}(vw)/(v). 

This is independent of a in the steady state and we 
shortly write (T cor (l|l)). In the same way as in the un- 
corrected case this yields 

oo 

<T(«l)> =T!=J2 ^(Tcor(0|l))(T col .(l|l)) ff 



(r cor (i|i)) 

l-<Tcor(l|l)> 



(54) 



Table IIVI compares simulation and mean field results 
of the mean occupation and the mean life time. They are 
all in very good agreement. 



TABLE IV. Simulation vs. mean field results for groups in 
the 2-cluster pattern for p — 0.025 







Si 


5 2 


S3 


mean occupation 


(n(v)) 3i 


0.993 


0.0004 


0.000 




UMFA 


0.993 


0.0003 


0.000 


mean life time 




6115 


0.017 


0.000 




TMFA 


6378 


0.014 


0.000 


occ. neighbours 


{n(dv)) Sl 


1.002 


10.94 


55.60 




«.(6Y>)mfa 


1.001 


10.94 


55.62 



VI. EVALUATION OF THE GENERAL CASE 

We notice that the solution of Eq. ([5]) can be very 
laborious, apart from very few simple special cases. In 



general we can determine the stable fixed points n* by 
iteration starting with some reasonable initial values n°. 

The suitable choice of the initial values is crucial, be- 
cause Eq. (J5J) often has multiple stable fixed points. How- 
ever, the basin of attraction is rather large for fixed points 
corresponding to dynamic patterns. 

Generic initial values, i.e. values near the stationary 
results of the respective simulation, converge to the sim- 
ulation results. 

For a nongeneric choice of n° the iteration may lead 
to results which do not correspond to the behaviour of 
the simulated system. For example, homogenous initial 
values n° g = n for all groups lead to a homogenous fixed 
point n* — n* . This is equivalent to an architecture with 
only one group. An analysis of one single polynomial 
equation like in subsection IIVBI can be possible in this 
case. [22I analyzes a similar situation. Symmetric initial 



conditions n° with n a „ — 



l d M +2-g 



lead to symmetric fixed 



points. This is due to the symmetry of the link matrix 

L g i = L dM+2 ~g,d M +2-l ) (55) 

which is inherited to the update map f(n). 

For generic initial conditions we computed the fixed 
points for various pattern modules and an interesting in- 

(2) 

terval of the parameter p on G 12 and [tL,tu] = [U 10] 
and compare the results with simulation data. 

In the simulations we obtained the mean occupation 

as 



1 



Ti -T 



E 



n t {v) 



(56) 



te(T ,Ti] 

the mean number of occupied neighbors of v 



n(dv) = V n t (dv), (57) 

J-i — Jo ' — ' 



t£(T ,Ti] 



and the mean life time 



t(v) 



1 



b(v) + n To 



(58) 



where b(v) is the number of births during the observation 
time, i.e. the number of new occupations of the node by 
the influx. Of course, b(v) + ut ^ must be fulfilled, 
otherwise t(v) has no meaning. 

In Fig. [8] we compare mean field results with the sim- 
ulation. In the left column we plotted isolines of his- 
tograms of n(i>), n(dv) and t(v) as a function of p, the 
simulations were started from empty base graphs. In the 
right column we used snapshots of the different patterns 
as initial conditions. 

In the simulations started from the scratch for differ- 
ent values of the influx p, the system evolved randomly 
and reached a steady state with a stable architecture. 
In the steady state we measured n(w), n(dv) and t(v) 
for all nodes during an observation time of 500, 000 it- 
erations. For each choice of p we created histograms 
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FIG. 8. (Color online) Statistical characteristics of nodes (mean occupation, mean life time, and mean number of occupied 
neighbors) vs influx probability p, in steps of Ap — 5/2 12 , obtained in mean field theory {solid lines) and by simulation. Left: 
Simulations with 500,000 iterations starting from the empty base graph. The thin gray lines are isolines (with increment 40) 
of histograms counting the frequency of nodes with a given characteristics after reaching the steady state, cf. Fig. 4 in [28| ]. 
Right: Simulations with 200,000 iterations started from a snapshot of the respective pattern. The symbols represent averages 
over nodes which are members of the same group. 

The mean field results for nodes of group S g in a pattern module of dimension oIm are given by solid lines labelled by (oIm; g)- 
They follow the ridges of the histograms (left) and coincide in most cases with the simulation of the prepared patterns (right). 
Groups with mean occupation near zero are not labelled. For further discussion see text. 
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of the frequency of nodes against one of the three sta- 
tistical node characteristics. In the figure we show the 
isolines of a compilation of all histograms in a range of 
p E [0,0.1]. The observed aggregations of nodes indicate 
the self-organization of the nodes into persisting groups 
as well as the statistical differences of the groups in de- 
pendence of p. 

In general, the mean field results show a good agree- 
ment with the simulation. For cLm — 2 we took the cor- 
relation correction into account, as described in Sec. IVl 

Mean field results for some patterns are not matched 
by a corresponding peak in the histograms. These pat- 
terns are seldom reached in simulations starting from the 
scratch, because their basin of attraction is too far away 
from the initial empty configuration or behind a high 
barrier which is difficult to be passed. 

The peaks which are not matched by a mean field re- 
sult may occur if a pattern does not persist over the ob- 
servation time. This affects the time average of the node 
statistics, see, for example, the left column at p « 0.0275. 

In the plots we also see that the character of the pat- 
terns changes with increasing p. For moderate influx, 
p 0.03, we identified almost static patterns associated 
with modules g?m = 2,4,6. 

For p » 0.03. . .0.08 there are no permanently occu- 
pied nodes, the pattern is dynamic, but a well defined 
group structure with du = H exists, cf. Fig. Q] This is 
the most interesting architecture as explained above, cf. 
also [HI, . The basin of attraction of the correspond- 
ing fixed point of Eq. [5] is large compared to the static 
patterns. 

For p gl 0.075 perturbations due to the random in- 
flux become so large, that the orientation of the pattern 
module in the base graph starts changing. Mathemati- 
cally, these reorientations are rotations [32|. After such 
a rotation there still exist 12 groups, but their identifi- 
cation by temporal averaging of node characteristics is 
impossible. The mean field theory does not consider ro- 
tations. It rather assumes that the nodes always remain 
in their groups. However, we can compare averages over 
the whole graph from the simulations with the mean field 
results averaged over all groups 

^ d M +l 

= ]g\ E I^K - (59) 

3=1 

where x g denotes one of the quantities n g ,r gi or n(dv g ). 
These averages are in good agreement, see plots of mean 
occupation and occupied neighbors in Fig. |8] 

We believe that these rotations of the pattern are the 
origin of the difference between mean field results and 
simulation in the mean life time plot at p « 0.85. Due to 
a rotation many previously occupied nodes become stable 
holes. Thus, significantly many nodes do not reach their 
usual life expectation. 

The fact that we do not observe all patterns in the 
simulations starting from an empty base graph, is a mo- 
tivation to prepare the initial conditions. In the right col- 



umn of Fig. [8] we compare the mean field results to node 
statistics from simulations with different predetermined 
initial architecture. The observation time was 200,000 
iterations. The picture proves that for the same value of 
p different patterns can persist, however, with different 
stability. The plotted lines of a pattern end at the value 
of p for which the pattern becomes unstable. 



VII. CONCLUSIONS 

We considered a minimalistic model to describe 
the random evolution of the idiotypic network which is, 
given very few model parameters, mainly controlled by 
the random influx of new idiotypes and the disappear- 
ance of not sufficiently stimulated idiotypes. Numerical 
simulations have shown that after a transient period a 
steady state is achieved. Depending on the influx and on 
other parameters, the emerging architecture can be very 
complex. Typically, groups of nodes can be distinguished 
with clearly distinct statistical properties. These groups 
are linked together in a characteristic way. 

We achieved a detailed analytical understanding of 
the building principles of these very complex structures 
emerging during the random evolution [281 ] . Modules of 
remarkable regularity serve as building blocks of the com- 
plex pattern. We can calculate for instance size and con- 
nectivity of the idiotype groups in perfect agreement with 
the empirical findings based on numerical simulations. 

In this paper we developed a modular mean field theory 
which allows to compute the statistical characteristics of 
a variety of patterns given their architecture. The re- 
sults are in very good agreement with the simulations. 
We thus have both structural and statistical informa- 
tion. The quantities usually considered in the context 
of networks, like degree distribution, cluster coefficient, 
centrality, etc., are coarser than those investigated here 
and they do not reveal the details of the architecture. 

For certain groups of nodes, namely singletons and the 
self-stabilizing core in static patterns, which decouple in 
good approximation from the rest, the mean field treat- 
ment can be simplified considerably. For the pattern 
consisting of pairs of occupied neighbors, the 2-cluster 
pattern, the naive mean field theory fails. Therefore, we 
extended the mean field theory to include correlations. 
Also for the cLm = H architecture correlations can be 
taken into account, which requires an intricate analysis 
of the building principles but is rewarded with an almost 
perfect match between theory and simulations [33| . 

The modular mean field theory can be applied to cal- 
culate statistical properties of other stationary networks 
with known modular architecture. 



Appendix: Stability 

Stability of patterns is a central and complex question. 
Here we restrict ourselves to discuss the stability of the 
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occupation state of a single node depending on the state 
of its neighborhood and the influx. We quantify this 
stability in terms of the probability that a node changes 
its occupation, i.e. to occupy empty nodes or to clear 
occupied nodes, respectively. 

In the following we calculate the probability to change 
the occupation of a node v given the number of occupied 
neighbors n(dv). we distinguish three cases. 

(1) n(v) = 0, n(dv) ^ tjj . This is an unstable hole, be- 
cause it can be occupied by the influx and will survive 
if n(dv) £ after the influx. The probability of 
occupation and survival is 

tu —n(dv) 

(k), (A.l) 

where p is the probability that v itself becomes occu- 
pied, and the sum gives the probability that the num- 
ber of occupied neighbors lies in the window, i.e. k £ 
[ti, — n{dv),tu — n(dv)] empty neighbors become occu- 
pied. PN.p(k) = ( I ^.)p k (l —p) N ~ k is the binomial distri- 
bution, where P/v,p(&) = if k < 0. 

(2) n(v) = 1, n(dv) < tu- This is an occupied node, in- 
cluding the case of an isolated node. The deletion of 
an occupied node requires that after the influx n(dv) is 
outside the window. For n(dv) < th the occupation of 
k £ [0, tjj — n(dv) — 1] or k £ [tjj — n(dv) + 1,k — n(dv)] 
empty neighbors is required. For n(dv) £ [tL,tu] the 
occupation of k £ [tu — n(dv) + 1, n — n(dv)] is needed. 
This yields 

(tL— n(dv) — 1 K-n(dv) \ 
E + E P,-n(dv)A k ) > 

k=0 k=tu-n(dv) + lj 

(A.2) 

where J2k=a • • • = if 6 < a. 

The changes of occupation in cases (i) and (ii) are 
possible within one iteration. There are more complex 
scenarios which require more iterations. All occupation 
changes which demand the deletion of neighbors need 
at least two iterations, for instance the occupation of a 
stable hole. Such occupation changes also involve the 
neighbors of the neighbors. As an additional sophisti- 
cation, we have to regard, that two neighbors of v have 
overlapping neighborhoods, they have at least the node 



v in common. 

(3) n(v)—0, n(dv)>tjj. This is a stable hole. During the 
first iteration the influx may increase the number of occu- 
pied neighbors by k £ [0, K — n(dv)], such that the num- 
ber of occupied neighbors becomes h(dv) — n(dv) + k. 
The probability for this process is P K - n (dv),p{k). In or- 
der to achieve n'{dv) ^ tjj, I nodes of dv must be emp- 
tied, I £ [h(dv) — tjj,n(dv)}. The probability to empty 
one occupied node w £ dv is given by P2(n(dw)), cf. 
Eq. (|A.2I) . However, the probability to empty I occu- 
pied nodes in dv is not simply n»=i ^ (n(dwi)), because 
the neighborhoods dwt are not disjoint. Neglecting this 
correlation, one can proceed with the second iteration as 
follows. We have after the first iteration n'(v) = and 
n'(dv) = n(dv) + k — I ^ tu- Thus, the second iteration is 
simply case (1), i.e. the occupation of an unstable hole. 
Collecting the probabilities of the consecutive steps for 
all possible choices of k and I gives then 

K-n(dv) 

P 3 (d(dv))= J2 P K-n{dv),p( k ) 
n(dv)-\-k 

x E Ell Mn(dw)) 

l=n(dv)+k-tv {C,}w£G'i 

x Pi(n{dv) +k — l), (A.3) 

where ^{cj ■ • ■ denotes the sum over all choices of I 
occupied nodes of dv which are emptied. 

For our parameter setting, and generally if <C tu <C 
k, we find Pi ^ Pi ^ P3. Stable holes are the nodes 
with the most stable occupation state. We expect that 
a pattern is stable, i.e. its nodes preserve their statisti- 
cal properties for a long time, if it has sufficiently many 
stable holes. 
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